CosMx_DBKero_BC_preprocessing_heatmap

Library import

library(SpatialExperiment)
library(data.table)
library(scater) # it imports scuttle
library(cowplot)
library(ggplot2)
theme_set(theme_bw())
library(matrixStats)
library(dplyr)
library(tidyr)
library(tibble)
library(arrow)
library(scales)
library(sf)
library(tmap)

Read-in function

readCosmxSPE <- function(dirname = dirname, 
                         countmatfpattern = "exprMat_file.csv", 
                         metadatafpattern = "metadata_file.csv", 
                         coord_names = c("CenterX_global_px",
                                         "CenterY_global_px")){
  
  countmat_file <- file.path(dirname, list.files(dirname, countmatfpattern))
  metadata_file <- file.path(dirname, list.files(dirname, metadatafpattern))
  
  # Read in 
  countmat <- data.table::fread(countmat_file)
  metadata <- data.table::fread(metadata_file)
  
  # Count matrix   
  counts <- merge(countmat, metadata[, c("fov", "cell_ID")])
  counts <- subset(counts, select = -c(fov, cell_ID))
  cell_codes <- rownames(counts)
  features <- colnames(counts)
  counts <- t(counts)
  
  rownames(counts) <- features
  colnames(counts) <- cell_codes
  
  # rowData (does not exist)
  
  # colData
  colData <- merge(metadata, countmat[, c("fov", "cell_ID")])
  
  spe <- SpatialExperiment::SpatialExperiment(
    assays = list(counts = counts),
    # rowData = rowData,
    colData = colData,
    spatialCoordsNames = coord_names
  )
  
  return(spe)
}

Setting directories

dirname <- "/NAS06/work/Spatial_omics/DBKero/CosMx_Breast/CosMX_data_Case2"
count_pattern <- "Run5810_Case2_exprMat_file.csv"
meta_pattern <- "Run5810_Case2_metadata_file.csv"
sample_name <- "CosMx_DBKero_BC"

CosMx dataset exploratory analysis for preprocessing

SpatialExperiment object creation and exploration

spe <- readCosmxSPE(dirname = dirname, countmatfpattern = count_pattern, metadatafpattern = meta_pattern)
spe <- scuttle::addPerCellQC(spe, subsets=list(NegProb = grep("^NegPrb", rownames(spe))))
names(colData(spe)@listData)

##  [1] "fov"                      "cell_ID"                 
##  [3] "Area"                     "AspectRatio"             
##  [5] "CenterX_local_px"         "CenterY_local_px"        
##  [7] "Width"                    "Height"                  
##  [9] "Mean.PanCK"               "Max.PanCK"               
## [11] "Mean.CD68"                "Max.CD68"                
## [13] "Mean.MembraneStain_B2M"   "Max.MembraneStain_B2M"   
## [15] "Mean.CD45"                "Max.CD45"                
## [17] "Mean.DAPI"                "Max.DAPI"                
## [19] "sample_id"                "sum"                     
## [21] "detected"                 "subsets_NegProb_sum"     
## [23] "subsets_NegProb_detected" "subsets_NegProb_percent" 
## [25] "total"

dim(spe)

## [1]  1010 59284

Selecting only numeric and integer variables from metadata then only variables of interest from the list

temp <- unlist(lapply(spe@colData@listData, function(x) is.numeric(x) | is.integer(x)))
num_variables <- c()
for (i in 1:length(temp)){
  num_variables[i] <- ifelse(temp[i] == TRUE, names(temp[i]), NA)
  num_variables <- num_variables[!is.na(num_variables)]
}
num_variables

##  [1] "fov"                      "cell_ID"                 
##  [3] "Area"                     "AspectRatio"             
##  [5] "CenterX_local_px"         "CenterY_local_px"        
##  [7] "Width"                    "Height"                  
##  [9] "Mean.PanCK"               "Max.PanCK"               
## [11] "Mean.CD68"                "Max.CD68"                
## [13] "Mean.MembraneStain_B2M"   "Max.MembraneStain_B2M"   
## [15] "Mean.CD45"                "Max.CD45"                
## [17] "Mean.DAPI"                "Max.DAPI"                
## [19] "sum"                      "detected"                
## [21] "subsets_NegProb_sum"      "subsets_NegProb_detected"
## [23] "subsets_NegProb_percent"  "total"

variables <- c("total", "detected")

Large-scale to small-scale

Approaching dataset preprocessing starting from global-scale, proceeding to local-scale, then fov-focused exploratory analysis

Sample map with fov grid

fov_pattern <- "Run5810_Case2_fov_positions_file.csv"
fovpos_file <- file.path(dirname, list.files(dirname, fov_pattern))

metadata <- data.table::fread(file.path(dirname, list.files(dirname, meta_pattern)))
fov_positions <- data.table::fread(fovpos_file, header = T)
fov_positions <- fov_positions[fov_positions$fov%in%unique(metadata$fov),]

# image theme
my_image_theme <- function(back.color="black",back.border=NA,title.col="white") {
  theme(panel.border = element_blank(),
        legend.key = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        panel.grid = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.title = element_text(color = title.col,hjust = 0.5,face = "bold"),
        plot.background = element_rect(fill = "transparent",colour = back.border))
}

fov_size <- 4256

ggplot() + 
  geom_point(data = metadata, mapping = aes(x = CenterX_global_px, y = CenterY_global_px), colour = "darkmagenta", fill="darkmagenta", size = 0.05, alpha = 0.2) +
  annotate("rect", xmin = fov_positions$x_global_px, 
           xmax = fov_positions$x_global_px + fov_size, 
           ymin = fov_positions$y_global_px, 
           ymax = fov_positions$y_global_px + fov_size,
           alpha = .2, color = "black",linewidth = 0.2) +
  geom_text(aes(x = fov_positions$x_global_px+fov_size/2,
                y = fov_positions$y_global_px+fov_size/2,
                label = fov_positions$fov), color="black", fontface = "bold") +
  ggtitle(sample_name) + 
  my_image_theme(back.color = "white", back.border = "white", title.col = "black")

Numeric/integer variable dependence on global position

Importing global coordinates from spatialCoords slot since they are not available in colData

colData(spe)$CenterX_global_px <- spatialCoords(spe)[,1]
colData(spe)$CenterY_global_px <- spatialCoords(spe)[,2]

Global coordinate scatterplot colored by variables showing spatial pattern

spe$CenterX_global_um <- spe$CenterX_global_px*0.18
spe$CenterY_global_um <- spe$CenterY_global_px*0.18  
spe$logtot <- log10(spe$total)

plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "total", point_size = 0.1) + ggtitle("Total counts per cell")

plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "logtot", point_size = 0.1) + ggtitle("Log-total counts per cell")

plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "detected", point_size = 0.1) + ggtitle("Detected features per cell")

spe$um_area <- spe$Area*0.0324
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "um_area", point_size = 0.1) + ggtitle("Area in μm per cell")

plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "AspectRatio", point_size = 0.1) + ggtitle("Width/Height ratio per cell")

spe$target_counts <- spe$total - spe$subsets_NegProb_sum
spe$target_detected <- spe$detected - spe$subsets_NegProb_detected
spe$tot_det_ratio <- spe$target_counts/spe$target_detected
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "tot_det_ratio", point_size = 0.1) + ggtitle("Target counts/detected features ratio per cell")

plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "subsets_NegProb_sum", point_size = 0.1) + ggtitle("Negative control probe counts/detected ratio per cell")

Check dependence on fov

Global coordinate scatterplot colored by binary color code: in black are the cells below the set variable threshold, in grey all the cells above

threshold <- 0.10
spe$var_color <- as.factor(ifelse(colData(spe)$total < quantile(colData(spe)$total, probs = threshold), "black", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Counts below 10% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.05
spe$var_color <- as.factor(ifelse(colData(spe)$total < quantile(colData(spe)$total, probs = threshold), "black", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Counts below 5% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

threshold <- 0.10
spe$var_color <- as.factor(ifelse(colData(spe)$detected < quantile(colData(spe)$detected, probs = threshold), "black", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Detected features below 10% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.05
spe$var_color <- as.factor(ifelse(colData(spe)$detected < quantile(colData(spe)$detected, probs = threshold), "black", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Detected features below 5% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

spe$area_um <- spe$Area*0.0324
threshold <- 0.10
spe$var_color <- as.factor(ifelse(colData(spe)$area_um < quantile(colData(spe)$area_um, probs = threshold), "black", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Area in μm below 10%") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.05
spe$var_color <- as.factor(ifelse(colData(spe)$area_um < quantile(colData(spe)$area_um, probs = threshold), "black", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Area in μm below 5%")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

threshold <- 0.90
spe$var_color <- as.factor(ifelse(colData(spe)$AspectRatio > quantile(colData(spe)$area_um, probs = threshold), "red", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(red = "red",lightgrey = "lightgrey")) + ggtitle("Area in μm above 90% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.95
spe$var_color <- as.factor(ifelse(colData(spe)$AspectRatio > quantile(colData(spe)$area_um, probs = threshold), "red", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(red = "red",lightgrey = "lightgrey")) + ggtitle("Area in μm above 95% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

threshold <- 0.10
spe$var_color <- as.factor(ifelse(colData(spe)$AspectRatio < quantile(colData(spe)$AspectRatio, probs = threshold), "black", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Width/height below 10% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.05
spe$var_color <- as.factor(ifelse(colData(spe)$AspectRatio < quantile(colData(spe)$AspectRatio, probs = threshold), "black", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Width/height below 5% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

threshold <- 0.90
spe$var_color <- as.factor(ifelse(colData(spe)$AspectRatio > quantile(colData(spe)$AspectRatio, probs = threshold), "red", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(red = "red",lightgrey = "lightgrey")) + ggtitle("Width/height above 90% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.95
spe$var_color <- as.factor(ifelse(colData(spe)$AspectRatio > quantile(colData(spe)$AspectRatio, probs = threshold), "red", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(red = "red",lightgrey = "lightgrey")) + ggtitle("Width/height above 95% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

threshold <- 0.10
spe$var_color <- as.factor(ifelse(colData(spe)$tot_det_ratio < quantile(colData(spe)$tot_det_ratio, probs = threshold, na.rm = T), "black", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Target counts/features below 10% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.05
spe$var_color <- as.factor(ifelse(colData(spe)$tot_det_ratio < quantile(colData(spe)$tot_det_ratio, probs = threshold, na.rm = T), "black", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(black = "black",lightgrey = "lightgrey")) + ggtitle("Target counts/features below 5% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

threshold <- 0.50
spe$var_color <- as.factor(ifelse(colData(spe)$subsets_NegProb_sum > quantile(colData(spe)$subsets_NegProb_sum, probs = threshold, na.rm = T), "red", "lightgrey"))
p1 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(red = "red",lightgrey = "lightgrey")) + ggtitle("Negative probe counts above 50% threshold") + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

threshold <- 0.75
spe$var_color <- as.factor(ifelse(colData(spe)$subsets_NegProb_sum > quantile(colData(spe)$subsets_NegProb_sum, probs = threshold, na.rm = T), "red", "lightgrey"))
p2 <- plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "var_color", 
            point_size = 0.1) + scale_color_manual(values = c(red = "red",lightgrey = "lightgrey")) + ggtitle("Negative probe counts above 75% threshold")  + theme(legend.position="none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1, p2, ncol = 2)

Heatmap with fovs colored by mean/median values of variable of interest

meta_df <- data.frame(colData(spe))
fov_cellcounts <- as.integer(table(colData(spe)$fov))
meta_df <- group_by(meta_df, fov) |> summarize(mean_total = mean(total), mean_detected = mean(detected), mean_areaum = mean(area_um), mean_aspectratio = mean(AspectRatio), mean_tot_det_ratio = mean(tot_det_ratio), mean_subsets_NegProb_sum = mean(subsets_NegProb_sum), median_total = median(total), median_detected = median(detected), median_areaum = median(area_um),median_aspectratio = median(AspectRatio), median_tot_det_ratio = median(tot_det_ratio), median_subsets_NegProb_sum = median(subsets_NegProb_sum))

meta_df <- cbind(meta_df, fov_cellcounts)

fov_positions <- left_join(fov_positions, meta_df, by = join_by(fov))

fov_size <- 4256
poly_df <- transmute(fov_positions, x_1 = x_global_px, y_1 = y_global_px, 
                     x_2 = x_global_px + fov_size, y_2 = y_global_px, 
                     x_3 = x_global_px + fov_size, y_3 = y_global_px + fov_size, 
                     x_4 = x_global_px, y_4 = y_global_px + fov_size)

str(poly_df)

## Classes 'data.table' and 'data.frame':   45 obs. of  8 variables:
##  $ x_1: num  151626 155892 155892 151626 137785 ...
##  $ y_1: num  13820 13820 18086 18086 23488 ...
##  $ x_2: num  155882 160148 160148 155882 142041 ...
##  $ y_2: num  13820 13820 18086 18086 23488 ...
##  $ x_3: num  155882 160148 160148 155882 142041 ...
##  $ y_3: num  18076 18076 22342 22342 27744 ...
##  $ x_4: num  151626 155892 155892 151626 137785 ...
##  $ y_4: num  18076 18076 22342 22342 27744 ...
##  - attr(*, ".internal.selfref")=<externalptr>

ls_list <- list()
for (k in 1:dim(fov_positions)[1]){
  ls_list[[k]] <- st_linestring(matrix(as.numeric(poly_df[k]), ncol = 2, byrow = T))
}

sf_data <- st_as_sf(fov_positions, coords = c("x_global_px", "y_global_px"), geometry = st_sfc(ls_list))

sf_data <- st_cast(sf_data, "POLYGON")

# plot(sf_data)
names(sf_data)

##  [1] "fov"                        "mean_total"                
##  [3] "mean_detected"              "mean_areaum"               
##  [5] "mean_aspectratio"           "mean_tot_det_ratio"        
##  [7] "mean_subsets_NegProb_sum"   "median_total"              
##  [9] "median_detected"            "median_areaum"             
## [11] "median_aspectratio"         "median_tot_det_ratio"      
## [13] "median_subsets_NegProb_sum" "fov_cellcounts"            
## [15] "geometry"

sf_data$fov <- as.factor(sf_data$fov)

heatmap_variables <- names(sf_data)[-c(1, length(names(sf_data)), length(names(sf_data))-1)]
for (var in 1:length(heatmap_variables)){
  heat_var <- heatmap_variables[var]
  p <- ggplot() +
  geom_sf(data = sf_data, aes(fill = .data[[heat_var]])) + 
  scale_fill_viridis_c() +
  theme(axis.title.x = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.title.y = element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid = element_blank()) +
  geom_text(data = fov_positions, aes(x = fov_positions$x_global_px+fov_size/2,
                y = fov_positions$y_global_px+fov_size/2,
                label = fov_positions$fov), color="white", size = 2) +
    ggtitle(paste0(heat_var))
  print(p)
}  

Boxplot of each of the numeric/integer variable distribution by fov, useful to see spatial patterns in data

variables <- c(variables, "logtot", "area_um", "AspectRatio", "tot_det_ratio", "subsets_NegProb_sum")

for (var in 1:length(variables)){
  n <- grep(paste0("^", variables[var], "$"), colnames(spe@colData))
  boxplot(spe@colData@listData[[n]] ~ spe$fov, ylab = colnames(spe@colData)[n], cex = 0.2)
  title(paste0(colnames(spe@colData)[n], "~fov"))
}

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 1 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 2 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 8 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 11 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 16 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 19 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 20 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 21 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 22 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 31 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 37 is not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 41 is not drawn

Per fov cell counts

fov_df <- data.frame("n_cells" = as.integer(table(spe@colData@listData$fov)), "fov" = as.factor(names(table(spe@colData@listData$fov))))
levels(fov_df$fov) <- order(levels(fov_df$fov), seq(1:length(fov_df$fov)))

ggplot(data = fov_df, aes(x = fov, y = n_cells)) +
  geom_col(fill = "#B0C4DE") + theme(legend.position ="none") +
  ggtitle("Number of cells per fov") 

Total/Area scatterplot

cell_df <- as.data.frame(spe@colData@listData)

ggplot(data = cell_df, aes(x = um_area, y = total)) +
  geom_point() + 
  ggtitle("Counts ~ μm area")

ggplot(data = cell_df, aes(x = um_area, y = AspectRatio)) +
  geom_point() + 
  ggtitle("Counts ~ μm area")

Detected/Area scatterplot

ggplot(data = cell_df, aes(x = um_area, y = detected)) +
  geom_point() + 
  ggtitle("Detected features ~ μm area")

AspectRatio/Area scatterplot

ggplot(data = cell_df, aes(x = um_area, y = AspectRatio)) +
  geom_point() + 
  ggtitle("AspectRatio ~ μm area")

Total/Area scatterplot faceted by fov

cell_df <- as.data.frame(spe@colData@listData)

ggplot(data = cell_df, aes(x = um_area, y = total)) +
  geom_point(size = 0.1) +
  facet_wrap(~fov) + 
  ggtitle("Total counts cell ~ μm area by fov")

Detected/Area scatterplot faceted by fov

ggplot(data = cell_df, aes(x = um_area, y = detected)) +
  geom_point(size = 0.1) +
  facet_wrap(~fov) + 
  ggtitle("Detected features ~ μm area by fov")

AspectRatio/Area scatterplot faceted by fov

ggplot(data = cell_df, aes(x = um_area, y = AspectRatio)) +
  geom_point(size = 0.1) +
  facet_wrap(~fov) + 
  ggtitle("AspectRatio ~ μm area by fov")